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We  have  been  investigating  the  application  of  unitary  quantum  lattice  algorithms  for  the 
evolution  of  quantum  turbulence  in  Bose  Einstein  condensates  (BEC)  [1]-[[20].  The  Hamiltonian 
for  the  cold  atom  BEC  trapped  in  a  magnetic  well  under  the  mean  field  approximation  is  given  by 


H  =  jd3x  r(-V2+V>  +  TjH 


(1) 


and  this  gives  rise  to  the  scalar  Gross-Pitaevskii  (GP)  equation  for  the  wave  function  (with 
h-1  =  m) 


.dii/  _  8H 


dt  dy/" 


-V2  +  Vext+g\y/\2 


¥ 


(2) 


The  superfluid  (mean)  velocity  of  the  BEC  can  be  discerned  from  the  continuity  equation 


^+V.(pvKC)=0  (3) 

where  the  density  p  =  yn//  * .  From  the  GP  equation,  we  thus  obtain 

Pvbec=/(^*v^-V^^*)  (4) 

From  the  Madelung  transformation  (for  winding  number  nw ) 

V  =  Jp  exp[inw(p]  (5) 

one  connects  the  phase  of  the  wave  function  with  the  mean  superfluid  velocity 

'sEC^Vcp  (6) 

While  this  leads  to  V  x  v6EC  =  0 ,  there  is  no  circulation  for  simply  connected  domains: 

j>cdl>\/BEC  =  0  :  if  domain  is  simplify  connected  (9) 


However,  it  has  been  found  experimentally  that  in  sufficiently  rotated  BECs  or  He4,  one  does  find 
quantized  circulations.  Hence  we  must  be  dealing  with  a  multi-connected  domain  so  that 


Fig.  1  Various  scalar  BEC  densities  under  rotation  (experiment)  :(a)  top  plot  -  non-rotating  BEC 
density;  (b)  plot-to-the-right :  rotating  BEC,  with  16  steady  state  quantized  vortices  (the  “holes"  are 
the  vortex  cores);  (c)  bottom-plot :  more  rapidly  rotating  BEC  with  70  vortices  in  lattice  formation, 
(d)  plot-to-the-left :  most  rapidly  rotating  BEC  with  130  vortices  in  lattice  formation. 
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i^-vc*0 


(10) 


with  topological  singularities 


p(*  c)  =  ° 


(11) 


and  quantization  of  circulation  about  these  singularities 
d)  dl»yacr=n  k 

J  BEC  w 


(12) 


where  k  is  the  quantization  of  circulation.  Thus,  for  rotating  scalar  BECs  one  finds  singular 
quantized  vortices  with  zero  density  at  the  cores.  This  is  what  is  seen  experimentally,  in  Fig.  1 . 

In  an  early  simulation  we  chose  as  an  initial  condition  a  winding  number  nw-  6  singular 

line  vortex  core.  This  is  a  highly  unstable  initial  condition  since  the  energy  of  the  6-fold  degenerate 
line  vortex  is  greater  than  6  times  the  energy  of  a  winding-number-1  line  vortex.  One  finds  vortex 
reconnection  and  tanglement  -  which  Feynman  defined  as  quantum  turbulence. 


Fig.  2  Tanglement  and  reconnection  of  a  degenerate  scalar  quantum  line  vortex. 

It  is  of  interest  to  discuss  the  relationship  between  quantum  turbulence  spectra  and  that  for 
classical  fluid  turbulence.  Typically,  classical  fluid  turbulence  is  discussed  for  incompressible 
turbulence,  with  the  Kolmogorov  /r5/3  kinetic  energy  spectrum.  Clearly,  since  compressibility 
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plays  an  important  role  in  quantum  turbulence  (since  a  scalar  singular  quantum  vortex  has  zero 
density  for  its  core  and  then  asymptotes  to  a  constant  density  in  a  healing  length)  there  should  be 
some  differences.  In  a  simulation  on  57603-grid,  we  found  the  total  kinetic  energy  spectrum  has  a 
triple  cascade 


Fig.  3  The  total  kinetic  energy  spectrum  for  scalar  quantum  turbulence 


In  this  (CAP-2  )  run,  we  did  not  have  enough  time  to  separate  the  incompressible  from  the 
compressible  kinetic  energies  (see  Fig.  5  for  such  a  separation  on  a  smaller  grid  of  30723)  -  but  the 
existence  of  a  triple  cascade  region  is  typical  of  multi-scale  physics.  For  example,  solar  wind  data 


fluid  scales  k  kinetic  scales 


Fig,  4  Magnetic  energy  spectrum  from  solar  wind  data  exhibiting  2  cascades:  fluid  (small  k)  scales 
exhibiting  the  k~ 5/3  spectrum,  and  a  kinetic  scales  spectrum  of  k-3  from  the  kinetic  (large  k)  scales. 


It  is  tempting  to  associate  the  small  k  k~m  -spectrum  in  Fig.  3  with  classical  turbulence  scales,  the 
intermediate  sharp  cascade  with  semi-classical  physics,  and  the  large  k  k-3  spectrum  with  quantum 
cascades.  It  should  also  be  pointed  out  the  significance  difference  between  classical  and  quantum 
turbulence:  in  classical  turbulence  it  is  essential  that  the  fluid  be  dissipative  -  i.e.,  non-zero 
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viscosity.  On  the  other  hand,  the  scalar  BEC  is  a  Hamiltonian  system  with  energy  conservation. 
Indeed,  using  the  Madelung  transformation  one  can  rewrite  the  GP  equation  into  the  more  familiar 
continuity  and  momentum  equations  and  then  identify  the  compressible,  incompressible  and 
quantum  kinetic  energies.  We  have  performed  detailed  spectral  analysis  on  a  quantum  turbulence 
run  on  30723-grid  for  both  non-degenerate  ( nw  - 1 )  and  degenerate  [nw- 2)  linear  vortex  cores. 


U)— I  4t000 _ k  (b) o=2  .1  fc  48000  k 


Fig.  5  The  energy  spectra  for  two  different  set  of  initial  conditions  on  the  linear  vortex  cores  :  (a) 
nw=  1 ,  and  (b)  nw-  2.  The  more  turbulence  state  (b)  gives  rise  to  ‘cleaner’  spectra.  The  blue 

curve  -  incompressible  kinetic  energy  spectrum,  red  curve  -  compressible  kinetic  energy 
spectrum,  gold  curve  -  quantum  kinetic  energy  spectrum 

Some  significant  outcomes  from  these  simulations: 

(a)  the  incompressible  energy  spectrum  Enc(/()«  ,or  *<200’ 

(b)  for  nearly  the  complete  spectral  range,  30  <  k  <  1000,  Enc(/()~  /T3.  There  is  basically  a 

single  cascade  in  the  incompressible  kinetic  energy  -  with  no  sign  of  a  Kolmogorov  region. 

(c)  both  the  compressible  and  quantum  kinetic  energies  exhibit  a  3-cascade  spectrum: 

small  k:  E  (k)~  k~513  ,  E  (k)~  k~513 

intermediate  k:  E  (k)~/T8  ,  E  (/c)~/T8 

large  k:  E  (k)~  k~3  ,  E  (k)~  k~3 

°  comp  \  /  ’  qu  \  ) 


Unitary  Quantum  Lattice  Gas  Algroithm  (QLG) 

Quantum  entanglement  is  the  backbone  of  quantum  computing  and  quantum  information 
theory,  with  qubits  as  the  building  blocks.  Unlike  the  classical  bit  (which  can  take  on  either  the 
value  ‘0’  or  ‘1  ’),  and  the  classical  Pbit  (which  can  take  on  a  probability  between  ‘0’  and  T),  the 
qubit  state  exists  on  the  unit  sphere  with 

|c3r)  =  r0|0)  +  r1|1) .  with  |70|2+|7i|2  =  1  (13) 
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where  y0  and  y1  are  complex  probability  amplitudes. 


Bit 

0 


Pbit 

o 


1  1  1 
Fig.  6  The  classical  bit,  the  probability  bit  and  the  qubit. 

For  quantum  entanglement,  one  requires  at  least  two  qubits: 

|c/i£/2)=roo|00)+roi|oi)+r10|io}+r11|ii).  04) 

We  then  apply  an  intertwined  set  of  non-commuting  unitary  collision  and  streaming  operators  so 
that  we  recover  the  finite  difference  form  of  the  scalar  GP  equation  (2)  on  a  simple  cubic  lattice.  In 
this  mesoscopic  approach  we  consider  the  spinor 

rYoi(x’t)^ 


a))  = 


7io(x>0 


(15) 


as  a  representation  of  the  scalar  wave  function  y/(x,t).  In  particular,  we  consider  the  unitary 
collision  operator,  based  on  the  3D  relativistic  Dirac  particle  dynamics  theory  of  Yepez, 

-  cose(x)  -/sine(x)  )  g{x)=E.la{x) 

w  4  24  1  1 


0D  = 


-/sin#(x)  cos#(x) 


(16) 


and  the  unitary  streaming  operator  (an  exponential  of  spin-1/2  matrics)  which  is  basically  just  a 
shift  operator  on  the  individual  spinor  components  of  Eq.  (15) 

. "  ' 


Ax,0 


))■ 


70i(x+Ax) 
^io(x) 


r10(x+Ax) 


(17) 


The  unitary  collision  operator  locally  entangles  the  qubit  amplitudes  at  each  nodal  site  while  the 
unitary  streaming  operator  spreads  this  entanglement  throughout  the  lattice.  The  non-commuting 
interleaved  sequence  of  collide-streams  give  rise  to  our  QLG  mesoscopic  finite  difference 
algorithm: 

|</>(x,f  +  Af)^  =  ()[x,o]|</>(x,f)^  (18) 

where  the  evolution  operator 

=  VIA W A  and  i0  =  S_,.„C0iMC0  (19) 


While  Eqs.  (16)-(20)  is  the  mesoscopis  algorithm  solved  on  the  supercomputers,  for  as  yet 
unspecified  q(x)  ,  it  will  form  a  representation  of  the  GP  Eq.  (2)  only  if  the  computational 
parameters  are  so  chosen  that  we  run  the  code  in  diffusion-like  ordering: 
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At~Ax2~e2  ,  and  Q— >e2Q  . 

Under  this  ordering,  the  spinor  time  advancement  yields: 

|^(x,f  +  Af)^  =  |^(x,f)^-/e2[-crxV2  +  Q]|^(x,f^  +  0(e4) 

To  move  to  the  GP  Eq.  (2),  we  take  the  zeroth  moment  of  the  spinor,  defining 
V '(x,f)=(  1  1  ).|  0(x,O)=ro1(x,f)  +  y1o(x,f) 
and  defining 

O(X’f)=Vrexf+9rk(X’0  f 

so  that 


/^  =  -VV  + 


V  ,  +  < 

ex/ 


¥  +  o(e2)- 


(20) 

(21) 

(22) 

(23) 

(24) 


There  are  some  important  remarks  that  must  be  made  about  the  QLG  algorithm: 

(a)  QLG  is  a  perturbative  algorithm  for  non-linear  problems, 

(b)  Since  it  is  a  lattice-based  algorithm  it  will  result  in  a  finite  difference  representation  of  the 
GP  Eq.  (24)  provided  the  parameters  are  so  chosen  to  yield  diffusion-like  ordering, 

(c)  It  might  first  be  tempting  to  consider  Trotter-like  decompositions  for  higher  order 

representations  of  the  collision  operator  CD  -  just  as  one  is  forced  to  do  in  Quantum  Monte 

Carlo  simulations  of  (equilibrium)  Hubbard  model  problems.  However  this  is  totally  counter 
productive  for  QLG.  In  QLG,  all  operator  effects  (from  both  streaming  and  further  collision 


steps)  are  required  to  give  higher  order  effects  over  the  simple  linear  expansion  for  CD . 
These  error  terms  are  2nd  order,  o(e2)  -  exactly  the  error  terms  arises  from  the  lattice 
discretization  effects  themselves. 

(d)  A  consequence  of  this  required  ordering  is  that  the  GP  order  parameter  must  itself 
be  small,  with  |y/(x,f|  =  0(e) . 

(e)  The  GP  Eq.  (24)  satisfies  both  conservation  of  particles  and  energy: 


jd2x\y/(x,tf  =N0= const.  ;  Jcf3x  y/*(-V2 +Vgxt)\i/  +  |g|^f 


=  E  -  const. 


(25) 


The  computational  QLG  algorithm  will  indeed  be  satisfying  the  diffusion-like  ordering  in  the 
Parameters  if  both  N0  and  E  are  conserved  during  the  run. 


It  is  straightforward  to  extend  the  QLG  to  spinor  BECs  or  coupled  BEC  systems:  the 
collide-stream  CD,S  interleaved  sequence  yields  a  second  order  accurate  lattice  representation  of 

the  operator  id/dt  +  V2 ,  while  the  potential  term  (in  the  collision  operator)  is  an  “inert”  term  giving 
the  nonlinear  couplings.  Our  original  QLG  algorithm  split  the  nonlinear  terms  from  the  collision 
operator:  in  this  algorithm,  the  collision  operator  was  simply  the  unitary  square-root-of-swap  gate 
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swap 


1-/  1+/ 
1+/  1-/ 


exp 


with  C2 

fy  ) 

*  01 

fy  1 
*  10 

swap 

*Toi  j 

(26) 


What  is  interesting  to  note  about  this  old  algorithm  is  that  C4  =  / ,  and  the  interleaving  non- 

&  Swap  ’ 


commutivity  of  J  =S  .  nC  .S  .  nC  .S  .  nC  .S  .  ,C  is  what  keeps  J \*l .  It  is  thus 

1  xO  -Ax,0  swap  +Ax,0  swap  -Ax,0  swap  +Ax,0  swap  r  xO 

somewhat  reminiscent  of  the  elements  of  a  Lie  algebra  which  lie  infinitesimally  close  to  the  identity 
operator. 


Benchmarking  Against  ID  Vector  Soliton  Collisions 

We  have  performed  some  benchmarking  runs  against  exactly  soluble  1 D  vector  soliton 
collisions.  Consider  a  coupled  set  of  ID  NLS  equations 


where  are  the  electric  field  polairzations  of  a  single  mode  propagating  down  an  optical 

birefringent  optical  medium,  with  cross  phase  birefringence  coefficient  B  that  couples  to  the  two 
polarizations.  For  these  so-called  Manikov  solitons  one  can  find  some  initial  vector  soliton 
amplitudes  that  will  yield  an  inelastic  collision:  i.e.,  for  a  unique  set  of  amplitudes  the  first  post¬ 
collision  amplitude  for  one  of  the  solitons  of  Q1  will  be  zero.  However  for  nearly  all  other  collisions 
there  are  always  2  solitons  for  each  polarization  (Figs  7  and  8). 


Fig.  7  The  1st  vector  soliton  collision.  Initially  both  polarizations  have  a  2-soliton  state  ( Q1  in  blue, 
Q2  in  red) :  one  a  left  traveling  soliton  (at  x  =  1000)  and  the  other  a  right  traveling  soliton  (at  x  = 
5000).  The  upper  plot  is  the  pre-collision  solitons  at  3  snapshots.  The  lower  plot  is  thelst-post 
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collision  soliton.  Note  that  during  this  1st  post-collision  state  there  is  destruction  of  the  left-traveling 
soliton  for  Q1 . 


This  is  possible  because  the  normalization  constraint 

J  c/x|Q.|2  =  const.  =  c.  ,  i- 1,2  (28) 

can  still  be  satisfied  by  a  2-soliton  for  each  polarization.  This,  of  course,  is  not  possible  for  a  scalar 
vector  soliton  with  its  1 -soliton  state  and  its  non-zero  normalization. 

In  Fig.  8  we  plot  the  maximum  amplitudes  of  the  2-solitons  for  each  of  the  polarizations  for  21 
collisions. 


max,|Oi(x.t)| 


Fig.  8  The  plot  of  the  local  maxima  for  the  2-solitons  for  each  polarization  ( Q1  in  blue,  Q2  in  red)  in 

21  collisions.  Note  that  only  for  the  1st  collision  is  the  left-traveling  soliton  annihilated  if  the  unique 
initial  set  of  amplitudes  appropriately  selected. 
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Parallelization  of  QLG  algorithms 

The  beauty  of  the  QLG  algorithm  will  not  only  run  on  a  quantum  computer  when  they 
become  available  because  of  its  unitarity  structure,  but  it  is  highly  parallelizeable  on  classical 
supercomputers.  We  comment  on  both  strong  and  weak  scaling  timings  :  in  strong  scaling  one 
fixes  the  grid  and  increases  the  number  of  cores.  For  ideal  parallelization,  the  wallclock  time  will 
decrease  by  a  factor  of  two  if  the  number  of  cores  are  doubled.  In  weak  scaling,  the  amount  of 
computation  is  kept  fixed:  the  factor  increase  in  the  lattice  grid  is  matched  by  the  factor  increase  in 
the  number  of  cores. 

On  the  Oak  Ridge  National  Laboratory  CRAY  XK7,  TITAN  ,  we  ran  up  to  288,000  cores  on 
grid.  For  strong  scaling  we  find  parallelization  efficiency  of  an  excellent  99.2% 


Cores 

Wallclock  (s) 

Ideal  (s) 

Speed-up  [ideal] 

150  000 

974~3 

974.?3 

1.00  [1.00] 

240  000 

613.45 

608.83 

1.59  [1.60] 

288  000 

511.24 

507.36 

1.91  [1.92] 

Table  1  Strong  scaling  on  TITAN,  Cl 

RAY  XK7  using  961 

JO3 

Weak  scaling  is  also  excellent,  as  we  scale  from  a  16003  grid  with  4096  cores  to  64003  grid  with 
262  144  cores. 


Grid 

Cores 

Wallclock  (s) 

Ideal  (s) 

16003 

4  096 

332.33 

332.33 

32003 

32  768 

334.08 

332.33 

48003 

110  592 

332.85 

332.33 

64003 

262  144 

334.48 

332.33 

Table  2  Weak  scaling  on  TITAN,  CRAYXK7 


The  problem  we  encountered  with  our  QLG  algorithm  on  the  CRAY  XK7  is  the  low  computational 
intensity  of  our  code.  Thus  we  could  not  use  the  GPUs  efficiently  on  the  CRAY  XK7  -  seeing 
saturation  effects  on  thousands  of  CPU/GPU  chips.  Basically,  as  originally  designed,  the  original 
GPU  chips  did  not  communicate  with  each  other. 


On  the  Argonne  National  Laboratory,  we  have  been  able  to  run  at  1.17  Peta  Flops  using 
2/3rds  of  the  full  IBM  BG/Q  ,  MIRA.  Strong  scaling  to  the  full  786  432  cores  is  excellent,  using 
96003 core 


#nodes 

Ranks  -  Mode  C32 

Time  (s) 

Speed-up  [ideal] 

16  384 

524  288 

817.1 

1.0  [1.0] 

32  768 

1  048  576 

389.7 

2.1  [2.0] 

49152 

1  572  864 

275.8 

3.0  [3.0] 

Table  3  Strong  scaling  on  IBM  BG/Q  MIRA  using  96003. 


The  weak  scaling  on  MIRA  is  again  excellent. 
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Grid 

Nodes 

Ranks -Mode  Cl  6 

Time  (s) 

8003 

32 

512 

323.4 

16003 

256 

4  096 

323.4 

32003 

2  048 

32  768 

323.7 

64003 

16  384 

262  144 

326.5 

Table  4  Weak  scaling  on  the  IBM  BG/Q  MIRA  using  1 MPI  rank/core 


Grid 

Nodes 

Ranks -Mode  Cl  6 

Time  (s) 

8003 

32 

2  048 

203.7 

16003 

256 

16  384 

197.1 

32003 

2  048 

131  072 

197.8 

64003 

16  384 

1  048  576 

209.6 

Table  5  Weak  scaling  on  the  IBM  BG/Q  MIRA  using  4  MPI  rank/core 


The  strong  scaling,  with  OpenMP  timings  on  a  51203-grid  on  32  (out  of  a  maximum  of  48)  racks, 
yields  a  parallel  efficiency  94.1%  with  a  1 .17  Peta  Flops: 


4  ranks 

8  ranks 

16  racks 

32  racks 

Wallclock  (s) 

406.11 

203.62 

106.58 

53.94 

Cores 

65  536 

131  072 

262  144 

524  288 

Parallel  Efficiency 

100% 

99.7% 

95.3% 

94.1% 

LI  d-cache 

88.64% 

89.13% 

89.11% 

88.79% 

DDR 

2.59% 

2.51% 

2.56% 

2.63% 

GFIops/node 

38.42 

38.35 

36.34 

36.12 

PFIops 

0.156 

0.311 

0.595 

1.174 

Table  6  Strong  scaling  on  the  IBM  Bi 

S/Q  MIRA  using  OpenMP  on  51 203 

Finally,  we  present  scaling  on  ERDC’s  SGI  ICE  X  TOPAZ  that  we  achieved  in  our 
successful  competition  to  be  given  CAP-2  allocation. 


Cores 

Wallclock 

Speed-up 

Ideal 

13  824 

3515.90 

1.00 

1.00 

27  000 

1792.63 

1.96 

1.95 

36  000 

1343.39 

2.62 

2.60 

52  920 

917.18 

3.83 

3.83 

74  088 

655.49 

5.36 

5.36 

Table  7  Strong  scaling  on  SGI  ICE  X  TOPAZ  (ERDC)  on  grid  42003. 
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Similarly,  for  weak  scaling  we  find  excellent  scaling. 


Grid 

Cores 

Wallclock  (s) 

Ideal  (s) 

6003 

216 

648.08 

648.08 

12003 

1  728 

647.48 

648.08 

18003 

5  832 

646.38 

648.08 

24003 

13  824 

647.27 

648.08 

36003 

46  656 

654.31 

648.08 

42003 

74  088 

655.49 

648.08 

Table  8  Weak  scaling  on  the  SGI  ICE  X  TOPAZ  (ERDC) 


QLG  for  SPINOR  BECs 

If  a  spin-2  87Rb  BEC  is  trapped  in  a  magnetic  well,  the  spin  degrees  of  freedom  are  frozen 
along  the  magnetic  field  and  the  T  =  0  mean  field  theory  of  this  BEC  is  given  by  the  scalar  GP  Eq. 
(24).  However,  if  the  spin-2  BEC  is  trapped  in  an  optical  well  then  the  spin  degrees  of  freedom 
require  a  spinor-GP  mean  field  theory  representation.  For  spin-2  this  results  in  5  coupled  GP  Eqs., 

with  the  spinor  wave  function 'F  =  [^_2  y/^  y/0  y/1  y/2  j  .  The  Hamiltonian  is  now  given 


by 


where  the  mean  density 

"M=XkM2-  (3°) 

m=- 2 

the  spin  vector  density 

fM=  X  n(x,f)u-(x.O  <31) 

m,n=- 2 

and  f  are  the  components  of  the  5x5  spin-2  matrices,  with  f  =(fx  ,fy  ,  fz  ),  and  A(x,t)  is 

m,n  r  r  ’  m,n  \  m,n  ’  m,n  ’  m,n )  ’  \  / 


the  singlet-pair  amplitude 

AM=XHVm(x,fV_m(x,f) 

m=-  2 


(32) 


When  there  is  a  phase  transition  from  normal  to  BEC  there  is  a  change  in  symmetry 
properties  under  parameter  variations.  It  was  finally  understood  that  the  spinor  wave  function  'F 
is  the  order  parameter  for  the  system.  Let  G  be  the  group  which  leaves  the  energy  functional  H, 
Eq.  (29),  invariant  and  H  the  invariant  subgroup  that  leaves  the  wave  function  'F .  The  factor 
group  G/H  is  the  order  parameter  manifold  with  topologically  invariant  defects  (“quantum 
vortices”). 
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The  physics  of  quasi  Nambu-Goldstone  modes  [Uchino  et.  a.  2010]  is  of  major  importance 
in  determining  the  low  energy  behavior  of  various  systems  with  spontaneous  symmetry  breaking. 
This  occurs  in  nearly  all  branches  of  physics,  particularly  with  phase  transitions.  The  signature  - 
i.e.,  the  so-called  “charge”  -  of  the  topological  defect/quantum  vortex  is  defined  by  how  the  order 
parameter  changes  in  a  closed  path  around  the  defect.  This  will  dictate  how  various  equivalence 
classes  of  vortices  will  interact  and  reconnect.  For  a  scalar  BEC,  the  order  parameter  manifold  is 
just  the  unitary  group  U(1),  the  circle  group  consisting  of  all  complex  numbers  with  absolute  value  1 
under  multiplication.  The  topologival  invariants  are  the  additive  group  of  integers  -  an  Abelian 

group.  The  quantum  vortices  are  singular  vortices  with  zero  density  at  their  cores:  .  p =|y/]2  =  0 , 

and  the  circulation  about  the  cores  are  quantized. 

However  for  BEC  in  an  optical  lattice,  with  the  order  parameter  being  the  spinor  wave 
function  much  more  complex  topological  defects  occur  -  in  particular  coreless  vortices  (skyrmions), 
vortices  with  rational  superfluid  circulation  quantization  (reminiscent  of  the  fractional  quantum  Hall 
effect),  and  even  unquantized  circulation. 

In  particular  we  have  followed  the  evolution  of  2  interlaced  skyrmions  in  an  external 
potential.  Here  we  will  be  dealing  with  non-singular  coreless  vortices  with  2  coupled  GP  equations 

for  the  order  parameter  wave  function  'P  =  |  y/+  y/  j  :  the  y/+-  core  (where  its  density  is 

zero)  is  a  vortex  ring  and  is  shown  in  green,  while  the  y/_  -  core  involves  a  standard  line  core 

region  which  then  folds  out  into  a  torus  (shown  in  orange).  We  report  on  the  simulation  of  3 
interlaced  skrymions,  Fig.  10.  In  Fig.  9  we  show  the  order  parameter  manifold  for  2  such 
skrymions 


Fig.  9  The  order  parameter  manifold  for  the  two  interlaced  skrymions,  with  the  component  y/+  in 
green  and  yr  in  orange. 

In  Fig.  10,  we  notice  that  the  time  evolution  of  the  y/+  -core  component  expands  from  a  vortex  ring 
into  a  capula  with  the  absorption  of  the  inner  structures. 
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Fig.  10  The  time  evolution  of  3  interlaced 
skyrmions  cores,  t//  (green)-component 
core  is  initially  a  vortex  ring  which  then 
evolves  into  a  single  coupla.  The  y/ 

(orange)-component  core  consists  of  a 
linear  line  core  which  folds  onto  itself 
toroidally.  As  time  evolves  the  inner 
skyrmions  are  absorbed  into  a  final  single 
skyrmion. 
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The  spectra  for  these  components  are  shown  in  Fig.  1 1 


lift 


|9.l  •  *  If  C 


(a)  spectra  for  vortex  rings,  (b)  spectra  for  linear  core  with  its  toroidal  envelope 

Fig.  1 1  the  change  in  the  spectrum  for  (a)  vortex  ring  component  and  (b)  linear  vortex  core  with  its 
toroidal  closure.  The  blue  curves  are  the  kinetic  energy  spectra,  while  the  red  curves  are  the 
quantum  energy  spectrum. 


•  Spin-1  BEC  VORTEX  RECONNECTIONS 

We  have  just  started  looking  into  vortex  reconnections  in  spinor  BECs.  First,  we  consider 
spin-1  BEC.  The  Hamiltonian  for  the  spin-1  BEC  is  basically  that  of  Eq.  (29) ,  except  there  is  no 
singlet-pair  interaction  (i.e.,  C2=0)  and  the  summation  over  m  now  runs  from  m  =  -1  to  +  1 . 

First  consider  the  case  where  there  is  no  spin-spin  interaction,  so  that  we  only  have  the 
standard  weak  Boson-Boson  coupling  interactions : 

ci=° 
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Consider  4  mutually  orthogonal  vortices  for  m  =  -1  , 0,  and  +1 .  We  find  that  if  there  are  no  spin- 
spin  interactions  among  the  m  -  components  there  are  no  vortex-vortex  reconnecetions:  in  each 
m-component  the  line  vortices  remain  constant  (but  with  a  slight  expansion  since  these  linear  line 
vortices  are  not  steady  state  solutions  of  the  coupled  spinor  GP  equations),  Fig.  12 


Fig.  12  The  evolution  of  4  parallel  vortices  in  the  y-direction  for  the  m=0  spinor  component.  (The 
m  =  1  spinor  component  has  similar  evolution  for  its  4  parallel  vortices  in  the  x  -  direction. 
Similarly  for  the  m  =  -1  spinor  component.  ) 


Now  suppose  for  the  initial  conditions  we  permit  only  for  the  m  =  -1  spinor  component  to  have 
within  it  orthogonal  vortices.  As  expected  one  now  finds  vortex-vortex  reconnection,  but  only  within 
this  spinor  component,  Fig.  13.  This  is  exactly  what  one  finds  for  scalar  BEC  vortex-vortex 
reconnection.  These  vortices  are  Abelian. 


Fig.  13  The  evolution  ofthem  =  -1  spinor  component  with  orthogonal  vortices.  One  now  finds 
vortex-vortex  reconnection.  There  is  no  reconnection  in  the  m  =  0  and  m  =  +1  spinor  components 
which  just  consist  of  4  parallel  vortices  ,  orthogonal  for  each  component. 


We  now  include  the  effect  of  the  spin-spin  interaction  term,  ci  i-  0  ,  in  the  Hamiltonian  Eq.  (29)  on 
the  vortex-vortex  reconnection  for  the  Abeiian  spin-1  BEC.  Consider  purely  orthogonal  4-vortices 
in  each  of  the  three  component,  as  in  Fig.  12  [in  Fig.  12  there  was  no  spin-spin  interaction  term  and 
no  vortex-vortex  interaction].  Initially  the  m  =  -1  four  line  vortices  are  parallel  to  the  x-axis  while  the 
m  =  0  four  line  vortices  are  parallel  to  the  y-axis  and  similarly  the  m  =  +1  four  line  vortices  along  the 
z-axis.  We  see  rung  vortices  form  in  the  reconnection  in  each  component. 
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m  =  -1 


m  =  0 


m  =  +1 


Fig.  14  The  effect  of  the  spin-spn  interaction  for  the  spin-1  BEC  with  initially  orthogonal  line 
vortices  -  this  is  the  counterpart  to  Fig.  12  in  which  there  was  no  spin-spin  interaction  in  the 
Flamiltonian  and  no  vortex-vortex  reconnection.  In  the  left  column  we  have  the  time  evolution  of 
the  reconnection  (  time  increasing  in  the  downward  vertically  direction).  There  is  phase 
information  on  the  vortex  lines :  blue:  zero  phase,  red:  phase  2*pi. 
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